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ABSTRACT 

We derive the density profile of cold dark matter halos using a self-similar accretion model. We show that 
if the dumpiness of the infalling matter is taken into account, then the inner density slope, 5 = dlogp/dlogr, 
is close to -L Compared with the density profiles predicted by different numerical simulations, we find that 
outside ^ 0. 1 % of the virial radius, our solutions agree best with the fitting formula proposed by Navarro et al. 
(2004), c/lnp/c/lnr = -2(r/r_2)", a ~ 0.17, where r_2 is a characteristic radius, inside which the density profile 
becomes shallower than isothermal (S < -2). 

Subject headings: cosmology: theory - dark matter - large-scale structure of Universe 



L INTRODUCTION 

In the standard ACDM paradigm dark matter halos evolve 
out of small fluctuations in the initial density field. Since the 
structure of virialized objects may depend on their formation 
history, many attempts have been made to link the observed 
density distribution of dark matter halos with the initial spec- 
trum of density fluctuation. High resolution N-body simula- 
tions indicate that density profiles of dark matter halos have a 
universal shape 

(cr/r„ir)(l+cr/ryiry 
whose dependence on the initial conditions is apparent only 
from the value of the concentration parameter c (Navarro, 
Frenk & White 1997; henceforth NEW). While the NEW 
model has been usually successful in reproducing the ob- 
served mass profiles of large galaxies and galaxy clusters, 
it seems to be inconsistent with the rotation curves of the 
dwarf galaxies. In addition, several other N-body sim- 
ulations (e.2. Moore et al. 1999; Fukushiae & Makino 2003; 
iRicotti 2003), while generally obtaining similar results at 
large radii, found different density slopes in the inner regions 
of halos, thus raising some doubts about the accuracy of the 
NEW model. 

The alternative approach to N-body simulations is 
provided by analytic modeling and in particular b 
self-simila r accret ion models (Eillmore & Goldreich 198 . 
iBertsching er 198^ . The choice of a particular set of ini 
tial conditions (the profile of an initial density perturbation 
is assumed to be a power law, AM/M oc r~^^) allows to re- 
cover an exact analytic solution for the final density distri- 
bution, that can be compared to the N-body results. Fur- 
ther, it can be shown that the appropriate value of e (and 
therefore the final density profile) is fixed by the spectrum 
of the initial d ensity fluctuat ions, e = (n + 3) /6, where n(k) = 
d\ogP{k)/dk jPeebles 19 801 Unfortunately, the very sim- 
plicity that is a strong point of the analytic models, can 
also lead to their downfall, if an important physical pro- 
cess has been left out. Purely radial infall that results from 
the evolution of a spherically symmetric density perturba- 
tion, produces an extremely steep inner density profile, 5 = 
d log p{r)/d log r < -2. As neither the observations no any 
of the N-body simulations produced a similarly steep den- 
sity distribution, attempts have been made to fix the self- 
similar models by giving collapsing halos some initial an- 



gular momentum ( White &Zaritskv 1991: Sikvie et al. 19971 
ISubramanian 2000: Nusser 200 1). Including angular momen- 
tum allows to reduce the steepness of the inner slope, down 
to -9e/(3e+ 1) = -(3n + 9)/(« + 5). However, in general the 
solutions, being dependent on the amount of added angular 
momentum, still did not reproduce the universal shape given 
by equation Q. 

In this paper we show that the discrepancy between the self- 
similar models and N-body simulations is naturally resolved 
by including the process of dynamical relaxation into the for- 
mer. We prove that even tiny deviations from spherical sym- 
metry are sufficient to make the dark matter particles velocity 
distribution isotropic in the inner region, provided the den- 
sity slope is steeper than -1. Furthermore, we show that for 
n < -I (e < 1/6) isotropic velocity distribution leads to the 
inner density slope S ^ -I. Since the spectral index, n{k), 
spans the range from -3 (for smallest halos) to -1 (for galaxy 
clusters), our results imply that if the dark matter gravity is 
the only force to be considered, then 6^-1 in all existing 
virialized halos. 

The paper is organized as follows. In §2, we show that 
density profiles, which are steeper than p oc r"', are inconsis- 
tent with unisotropic velocities at the central region. In §3, 
we derive a density profile for a self-similar accretion model 
that includes dynamical relaxation process. In §4, we com- 
pare our results with the N-body codes. We summarize our 
conclusions in §5. 

2. DYNAMICAL RELAXATION 

In the standard ACDM paradigm the structure formation is 
hierarchical. The smallest halos that form first subsequently 
merge into larger and larger objects. If dynamical friction be- 
tween merging dark matter clumps is weak, dumpiness and 
radial orbits can be retained for a long time. On the other 
hand, if dynamical friction is strong, the clumps are destroyed 
and the velocities of dark matter particles become isotropic. 
As we show below, the latter is always the case at the cen- 
tral regions of halos with sufficiently steep (6 > -I) density 
profiles. 

Consider a halo of mass M/, and radius /?/, formed by 
smaller clumps with the typical mass md- The clumps char- 
acteristic relaxation time is 



trei ~ (TmciG^mlVd'^log[Mh/m,i])~ 



(2) 



where rid ~ p/nid and Vd are the clump number density and 
characteristic velocity. The kinetic energy of clumps must be 



2 



less than their potential energy at the center, Vj/2 < GMi,/Rh- 
This sets an upper limit on the relaxation time 



become 



trel < 



ml 



■KHcl 



•tl,R'Jhog[MH/mc,] 



(3) 



Since the radial infall velocity of a clump is less than ~ 
yj2GMh/Rh, the time required for a clump to reach the cen- 
ter from a distance r, f^r, exceeds r/ yJlGMijRj,. For radial 
orbits to be sustained, the crossing time must be less than the 
relaxation time 



K — < 



ml 



iTmci\og{Mh/mci\R 



[pr]- 



(4) 



Obviously, if the density profile is steeper than p oc r~^, then 
close to the center (i.e. for small values of r) the above in- 
equality would be violated, implying that in such a case dy- 
namical friction must be strong and radial orbits cannot be 
sustained. Therefore, in the case where purely radial infall 
leads to a density slope below —1, dynamical relaxation must 
either isotropize the velocities or flatten the central density 
cusp to p (X r~'. 

3. SELF-SIMILAR INFALL 

We now proceed to calculating the actual density distri- 
bution of the collapsing halos, including the dynamical re- 
laxation process. In the Einstein-de Sitter model (which is 
an excellent fit at high redshifts) a system evolving from a 
power law density perturbation (AM/M oc r~^'^) over a time 
approaches self-similarity, which allows to use dimensionless 
variables for mass profile, density, bulk velocity and energy 
density of dark matter particles. 



m(r,0=M(A)^, 


(5) 




(6) 


v(r,t) = V(X)'f, 


(7) 




(8) 



where A = r/r,,,> and f is a Hubble time. The virial radius, 
r,.„ , which is customarily defined as a radius inside which the 
density contrast is IStt^ « 178, grows as f'', where r] = 2(3e + 
l)/9e. 

At a fixed A the relaxation time (equation Q) scales as f, 
which implies that the relaxation process does not break self- 
similarity. However, for the sake of simplicity, rather than 
evaluate the impact of dynamical relaxation on particles ve- 
locities at each point, we assume that at large radii the infall 
is purely radial, neglecting small tangential motions. Close 
to the center we assume that the velocities are completely 
isotropic. Furthermore, we assume that an abrupt transition 
between the two regimes occurs at some point. This transition 
is analogous to the shockheating of the infalling gas (though 
for collisionless dark matter the width of the transition region 
can be significantly larger than the width of the shock, which 
roughly equals baryon mean free path). This approach allows 
us to treat the dark matter as a dissipationless fluid whose dy- 
namics is determined by continuity, Euler and adiabatic equa- 
tions. Cast into the dimensionless variables, these equations 



iV-7jX)D'+ (^ + V'-2^D = 0, 
W 2M 

{TJ-l)V + (V-7j\W'=-—-^, 

5D' U'\ , 2 

(y-7?A) = 2rj--, 

3D n J ' ' ' 3 ' 

M' = 3X^D, 



(9) 

(10) 

(11) 
(12) 



If all non-adiabatic cooling and heating processes can be ne- 
glected, the above equations are equally applicable to the 
baryonic and dark matter fluids (Chuzhov & Nusser 2000). 

Outside the shock the radial motion of the infalling mat- 
ter can be traced analytically ( Bertschinaer 1985). The jump 
conditions at the adiabatic shock are 



V^ = T]Xs+- 

D^ = 4D~, 
P^=^D-(V-r,X,,f 



(13) 
(14) 
(15) 
(16) 

where the superscript minus and plus signs refer to pre- and 
post- shock quantities. The location of the shock. A,, which is 
initially unknown, have to be found by integrating equations 
(I9l- (I12> inwards. Too small values of A, produce a non-zero 
mass at A = 0, while too large values of Aj give a singularity 
at some A > 0. A unique value of A, can therefore be recov- 
ered by requiring the mass and the infall velocity to be zero 
at A = 0. Typically, we find that the value of Aj is close to 
the location of the first dark matter shell-crossing, which is 
obtained when dynamical relaxation is neglected. 

The asymptotic analysis of equations (I9t-( ll2t reveals the 
existence of two different types of solutions. For -3 < n < -2 
(e < 1 /6) the matter energy density is finite everywhere and 
close to the center the density slope goes to 5 = -3(n + ?)/(« + 
17), thus spanning the range between -6/7 « -0.86 and -1. 
Since completely isotropic velocities, which we assumed at 
the center, require 5 < -1, it is unclear whether S > -I can 
actually be obtained. However, the convergence of S to its 
asymptotic limit at the center is rather slow (at r = 0.01r,,,> 
the slope is still ^ -1.3 and -1.2 for n = -2 and -2.5, respec- 
tively), (5 + 3(n + 7)/(n+ 17) cx 1/logr. Thus, except for ex- 
tremely small values of r, our solutions should remain valid. 

For n > -2 the matter energy density goes to infinity at the 
center and the density slope goes to -(3n + 9)/(n + 5).' How- 
ever, since n > -2 is obtained only on scales, whose evolution 
is still in the linear regime, these solutions do not apply to any 
existing halos. 

4. SELF-SIMILAR SOLUTIONS VS N-BODY CODES 

We find that for 0.01 < r/ r^r < 1 our solutions can be well 
approximated by the NEW fit (Figures ^ and |2j. The dis- 
crepancy of order 10% in the circular velocity profile, which 
is seen at the outer part, where "shock" takes place, is quite 
expected, given our simplistic modeling of this transitional 
region. Closer to the center ( 0.0001 < r/r,,„ < 0.01), our so- 
lutions predict steeper density slope than the NEW fit (though 
not as steep as predicted by Moore et al. (1999)). This is 

' Accidentally, when the dynamical relaxation process is neglected, the 
same slope is also obtained for n> \. 
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Fig. 1 . — Circular velocity [Vc = \/ Gm(r)/r] profile for (i = — 2 (e = 1/6). 
The dotted and the solild lines show the curves, respectively, with and without 
dynamical relaxation included. The dashed line shows the NFW profile for 
c = 4.5. 
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Fig. 2.— Circular velocity [Vc = Gm{r) / r] profile for n =-2.5 (e= 1/12). 
The dotted and the solild lines show the curves, respectively, with and without 
dynamical relaxation included. The dashed line shows the NFW profile for 
c = 5.3. 



consistent with the later simulation results of Navarro et al. 
(2004), who suggested an improved fitting formula for the 
density profile 



In 



p(f) 

P-2 



— 1 -1 

r-2 



(17) 



where r_2 is the radius inside which the density slope is above 
-2 and p_2 = p{r-2)- We find that at small radii the equation 
in\ is a better approximation to our results (see Figure |3} 



and the values of a that produce the best fits to our solutions 
in the range 0.001 < r/r_2 < 1 (0.16, 0.21 and 0.22, respec- 
tively, for n = -2, -2.5 and -2.75) are close to a = 0.17 found 
by Navarro et al. (2004). However, at r < 0.001r_2 (which 
has not been resolved by the N-body codes) the agreement 
between the equation J17> and our results breaks down. 

5. SUMMARY 

We have shown that the results of high-resolution N-body 
simulations of dark matter infall can be quite accurately re- 
produced by a simple self-similar accretion model. Moreover, 
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Fig. 3. — The dark matter density versus radius. The dotted, dashed- 
dotted and dashed lines correspond, respectively, to n = — 2, -2.5 and —2.75. 
The solid line corresponds to the fitting function (Eg. 1171 with a = 0.17. 



the self-similar accretion model, being unlimited by compu- 
tational constraints, may have an advantage over the N-body 
codes in reconstructing the density profile in the central region 
of dark matter halos. Naturally, it should be kept in mind that 
close to the center of a halo, a stellar disk or a supermassive 
black hole may strongly affect the density profile, making in- 
adequate any model or code that rely exclusively on the dark 
matter gravity. 

It is interesting to note that the two threshold values of n, 
which separate different types of self-similar solutions, have a 
special significance in the real Universe. Thus n = -2, which 
separates the solutions with finite and infinite central energy 
density, is also the largest possible value of n on scales of 
existing virialized halos. The « = 1 threshold, which separates 
the solutions with finite and infinite gravitational potential in 
the center, seems als o to be the largest value of « on all scales 
jSnergeletal.20n(i). 

LC thanks the McDonald Observatory for the W.J. McDon- 
ald Fellowship. 
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